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The quasi-bound electromagnetic modes for the arrays of nanoholes perforated in thin gold film 
are analyzed both numerically by the rigorous coupled wave analysis (RCWA) method and semi- 
analytically by the coupled mode method. It is shown that when the size of the nanohole occupies 
large portion of the unit cell, the surface plasmon polaritons (SPPs) at both sides of the film are 
combined by the higher order waveguide modes of the holes to produce multipole surface plasmons-. 
coupled surface plasmon modes with multipole texture on the electric field distributions. Further, 
it is revealed that the multipole texture either enhances or suppresses the couplings between SPPs 
depending on their diffraction orders and also causes band inversion and reconstruction in the 
coupled SPP band structure. Due to the multipole nature of the quasi-bound modes, multiple dark 
modes coexist to produce variety of Fano resonance structures on the transmission and reflection 
spectra. 

PACS numbers: 42.79.Gn, 73.20.Mf, 42.25.Bs, 78.66.Bz 


I. INTRODUCTION 

Since the discovery of the extraordinary optical trans¬ 
mission (EOT) of sub-wavelength hole array systems by 
Ebbessen [1], plenty of researches have been devoted to 
elucidate this mechanism (for recent review, see [2]). To 
date, the existing theories reveal that this phenomenon 
comes as a result of resonances with quasi-bound elec¬ 
tromagnetic (EM) modes localized around the metal film 
[2, 3]. Although sub-wavelength holes transmit light with 
a very low efficiency in general [4, 5], the energy of the 
incident light is carried efficiently to the opposite side 
of the film by exciting a quasi-bound mode in the film. 
Thus the EOT can be considered as an effect of resonant 
tunneling of light via quasi-bound modes. 

The quasi-bound modes in metallic nanohole arrays are 
regarded as coupled surface plasmon modes, in which the 
surface plasmon polaritons (SPPs) at both sides of the 
metal film are combined by the waveguide modes in the 
nanoholes. The contributions of the nanohole array are 
two folds. Eirstly, the SPP at each side of the film is mod- 
ihed due to the presence of the holes. This means that the 
periodic hole array does not only open band gaps in the 
SPP band structure but also produce confinement of the 
EM fields inside the hole array region. Indeed, it is known 
that geometrically induced surface EM modes, also called 
as “spoof surface plasmons”, appear in perforated PECs 
whose flat surfaces do not support SPPs [6]. Secondly, 
the (spoof) surface plasmons at both sides of the film 
are coupled through the evanescent fields of the waveg¬ 
uide modes in the holes and form two separate “plasmon 
molecule” levels [7]. In the symmetric environment of 
equal dielectric constants in the regions of incidence and 
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transmission, the two molecule levels correspond to the 
bonding and anti-bonding modes to which in the (anti- 
) bonding mode the SPPs at both sides of the film are 
combined (anti-) symmetrically and has a (higher) lower 
energy. 

The propagation constant, q^z, of the fundamental 
waveguide mode of the nanohole (HEn mode in the case 
of cylindrical hole) is the key parameter to determine 
the property of the coupled SPP modes. Such that, if 
goz becomes nearly zero at the cutoff frequency, then the 
confinement of the EM field becomes so large that the dis¬ 
persions of the coupled SPP modes flatten. Also, when 
Qoz — 0) then the condition corresponds to the zero or¬ 
der Fabry-Perot resonance and the resulting transmission 
maxima are much more like cavity resonance of a single 
hole [8, 9]. Therefore, the coupled SPPs around the cut¬ 
off frequency may be considered as cavity arrays that 
are weakly coupled by the SPPs, instead. On the other 
hand, when the imaginary part of qoz becomes large, the 
coupling between the waveguide mode and the SPPs be¬ 
comes weak and the coupled SPPs are almost decoupled 
and behave like (spoof) surface plasmons on two discon¬ 
nected metal surfaces. In this case, the EOT has the 
same origin as Wood’s anomaly [10, 11] which is deemed 
as a result of the EM waves propagating along the metal 
surfaces, such as the waves diffracted parallel to the sur¬ 
face of a PEC (Rayleigh anomaly [12]) or by the SPPs 
propagating on a real metal surfaces (plasmon anomaly 

[13]). 

However, there is another important ingredient of the 
EOT that has not been considered much so far. When 
the size of the hole occupies large portion of the unit 
cell, the higher-order waveguide modes would play an 
important role. These higher-order waveguide modes as 
shown in Figure 1 have multipole nature that weakens 
the coupling with incident light that has a dipole nature 
and hence leads to having a minor contribution to the 
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Figure 1. (Color online) Electric field distributions for the 
waveguide modes of a nanohole. 


transmission spectra. Nonetheless, it is known that the 
combination of a bright dipole mode and a dark multipole 
mode yields to a sharp peak structure in the transmission 
spectra via Fano resonance that focuses much attention 
because of the expectation of potential applications such 
as the high-resolution bio-sensors [14, 15]. Therefore, it 
is expected that a variety of sharp resonance structures 
will appear in the transmission and reflection spectra due 
to the existence of higher-order waveguide modes [16]. 

Moreover, it is expected that the multipole modes in 
the nanoholes are combined with each other by the SPPs 
on the film surfaces. In other words, the SPPs are com¬ 
bined by the multipole modes to produce novel type of 
hybridized bound modes which may be called as mul¬ 
tipole surface plasmons (MSPs). The band structure of 
MSPs should have a distinct feature that is rarely seen in 
the lattice of dipoles such as the natural crystal of atoms 
with higher order of atomic orbitals that generates vari¬ 
ous structures of diverse physical properties. 

In this paper, we investigate the contribution of mul¬ 
tipolar waveguide modes of nanoholes to the creation of 
quasi-bound EM modes for the arrays of nanoholes per¬ 
forated in thin gold film. In order to confirm the exis¬ 
tence of MSPs, we have performed a thorough numeri¬ 
cal simulation by using the rigorous coupled wave analy¬ 
sis (RCWA) method, which gives almost exact results of 
transmission, reflection and absorption spectra and also 
dispersion of quasi-bound modes. The contribution of 
the higher-order waveguide modes and their multipolar 
field distribution on the band structure of MSPs are fur¬ 
ther analyzed through the use of the coupled mode (CM) 
method. 


II. METHOD 

We adopt two different theoretical approaches, namely, 
the RCWA method ]I7, 18] and the CM method [2, 19]. 
The accuracy of the results obtained from the use of our 
homemade numerical codes are checked through compar¬ 
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Figure 2. (Color online) Schematic definition of the coefficient 
vectors for the EM modes entering and leaving the interfaces 
of layers. 


ison with those obtained by the finite-difference time- 
domain (FDTD) method [20] with the help of a freely 
available software package ]21]. In this section we present 
the basic ingredients of these two methods and explain 
their efficiency in obtaining the dispersion relations of the 
quasi-bound modes. 


A. RCWA method 

The RCWA method is applicable to systems that can 
be separated effectively into several planar layers. Each 
layer is invariant under translation in the thickness direc¬ 
tion which is assumed to be in the z-direction and is peri¬ 
odic in two other spatial directions {x- and y-directions). 
The fields in each layer are expanded in Fourier-Bloch 
modes that propagate or decay either forward or back¬ 
ward in the z-direction. We can derive the solution of a 
stack of layers in a scattering matrix formalism. Scatter¬ 
ing matrix S provides a relation between the coefficient 
vectors Ai and B 2 for the input modes and the coeffi¬ 
cient vectors A 2 and Bi for the output modes as shown 
below: 



( 1 ) 


where the vectors Ai (Bi) with i being 1 or 2 contain 
a set of coefficients that define the amplitude and phase 
of the EM modes propagating in the -j-z (—z) directions 
at the top (i = 1) or bottom (i = 2) interfaces (see 
Fig. 2). In this method, all the EM modes are expressed 
by the coefficient vectors with the same number of di¬ 
mensions A^fb that are determined by the number of 
Fourier-Bloch modes used in the expansion. Likewise, 
the scattering matrix S can be constructed iteratively for 
the stack of any number of layers considering the propa¬ 
gation in each layer and the boundary conditions at the 
interfaces [18]. Conversely, the truncation of Fourier se¬ 
ries with a finite value of A^fb causes serious problem of 
poor convergence when there is any discontinuous jump 
in the permittivity distribution. Although the conver¬ 
gence can be improved using the so-called Fourier factor¬ 
ization rules [17], the improvement is not enough in the 
case of metallo-dielectric structures since the difference 
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between the permittivity of the metallic and the dielec¬ 
tric area is quite large. Therefore, we adopt the scheme of 
matched coordinates developed by Thomas Weiss [18], in 
which the spatial resolution is increased adaptively close 
to the jump discontinuities with the application of the 
Fourier factorization rules by introducing a coordinate 
transformation technique developed at the beginning of 
transformation optics [22, 23]. (Figure SI in the Supple¬ 
mental Material [24] shows the coordinate we used for 
the triangular lattice of nano holes). 

From the comparison with the results obtained by the 
FDTD method, we see that the RCWA method gives 
almost exact results for the transmission and reflection 
spectra (see Figure S4 in the Supplemental Material). 
Throughout our calculations, we used TVfb = 29 x 29 
Fourier-Bloch modes. Figure S5 in the Supplemental Ma¬ 
terial shows that converging behavior is obtained around 
this number of Fourier-Bloch modes, and the expected 
error is small enough. The computational time required 
in the RCWA method is significantly shorter than that 
required to obtain quantitative results by the FDTD 
method. 


B. CM method 


express the orthogonality condition for the general hybrid 
modes having both TE and TM components. 

Based from the results of the comparison made with 
the RCWA and FDTD methods, we have shown that our 
modifications have improved the CM method enough to 
give the quantitative prediction of the whole resonant 
structure in the transmission and reflection spectra with 
a computational time that is four orders of magnitude 
shorter than that required by the FDTD method (see Fig¬ 
ure S2-S4 in the Supplemental Material [24]). Through¬ 
out our calculations, we used 91 Fourier-Bloch modes 
whose wavenumbers are not exceeding 5 times the mag¬ 
nitude of the basic reciprocal lattice vector. Figure S6 
in the Supplemental Material shows that converging be¬ 
havior is obtained around this number of Fourier-Bloch 
modes, and the expected error is small enough to obtain 
quantitative results. 

By the mean field approximation for the admittance 
operator, we can derive a similar expression as the origi¬ 
nal CM equations: 
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In the CM method, the EM fields in the metallic 
nanohole array region are expressed only by the waveg¬ 
uide modes of the nanoholes. Imposing the continuity 
condition of EM fields at the openings of the holes and 
the surface impedance boundary conditions (SIBCs) [25] 
at the surfaces of the metal film, we can derive a coupled 
system of equations for the coefficient vectors as follows: 
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where A (B) is the coefficient vector for the waveguide 
modes propagating in the +z (—z) direction at the top 
(bottom) interface. The matrix G^ controls the EM cou¬ 
pling between waveguide modes at the interfaces, and Qz 
denotes the diagonal matrix whose each diagonal element 
qaz is the propagation constant associated with mode a. 
The thickness of the metal film is denoted by h while the 
vector I takes into account the direct initial illumination 
over the waveguide modes. The mathematical expres¬ 
sions for the different magnitudes can be found in the 
Appendix. 

At this point, we extend the original formulation [2] in 
two aspects in order to improve its quantitative aspect: 
(i) waveguide modes for a real metal with finite permit¬ 
tivity are used for the calculation of the overlaps of modal 
wavefunctions with plane waves, while PEC waveguide 
modes were approximately used in the original treatment 
and (ii) the admittance operator is introduced in order to 


where E {E’) approximately represents the modal am¬ 
plitude of the electric field at the top (bottom) interface. 
The mean field approximation deviates considerably from 
the exact calculation in the metal region near the hole 
edge as we can see from the in-plane magnetic field dis¬ 
tributions shown in the Figure S7 of the Supplemental 
Material. However, the overlaps with plane waves are 
predominantly determined by the integral in the hole re¬ 
gion, the influence on the spectra is not significant as we 
can see in the Figure S3 of the Supplemental Material. 


C. Quasi-bound mode 

The occurrence of resonant transmission implies the 
existence of a quasi-bound mode in the scattering object. 
Quasi-bound modes are solutions of Maxwell’s equation 
that can oscillate in time and hold EM energy within the 
object for a considerable period of time in the absence 
of incident light. Formally, we would have to derive a 
non-vanishing output for zero input. 

In the RCWA method, this means that the quasi¬ 
bound modes would correspond to the output vectors 
A 2 and Bi of Eq. (1) for zero input vectors Ai = 0 and 
B 2 = 0. In other words, the vectors Ai and B 2 must be 
enlarged infinitely by the scattering matrix S. Numeri¬ 
cally, the search for quasi-bound modes can be fulfilled by 
seeking the infinite singular values of S using the routine 
of singular value decomposition. 

While in the CM method, the quasi-bound modes 
would agree to the non-trivial solutions of Eq. (2) for 
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Figure 3. (Color online) Schematic diagram of a triangular 
lattice of nanoholes perforated in a thin gold film. 

zero input vector 1 = 0, and thus we seek the eigen¬ 
vectors of zero eigenvalue (i.e. null space) of the matrix 

G. 

In searching for the quasi-bound modes, we must take 
special care of the boundary condition at infinity. Quasi¬ 
bound modes are not true bound modes since they have 
finite lifetime as they eventually lose energy due to the 
weak emission in the absence of incident light. There¬ 
fore, we must seek the EM eigenmodes in the lower half 
of the complex plane of angular frequency cu, adopt¬ 
ing the out-going-wave boundary condition. This means 
that one of the branches of the multi-valued function kz 
must be chosen obeying the condition that Refcz > 0 for 
IRefc^l > llmfc^l or Imfc^ > 0 for IRefc^j < jlmfe^j. Cer¬ 
tainly, this condition becomes nonphysical such that (i) 
mode fields explode at infinity (i.e. Imfcz < 0) for the 
modes predominantly propagating to infinity and (ii) en¬ 
ergy flows come from infinity (i.e. Kekz < 0) for the 
predominantly evanescent modes [26]. 


III. RESULTS 

A. Transmission, Reflection and Absorption 
Spectra 

The system that we are concerned in this paper is a 
triangular lattice of nanoholes perforated in a thin gold 
film with a thickness oi h = lOOnm as shown in Figure 3. 
The film is surrounded with water. The diameter of the 
hole is assumed to be 300nm and the period of array, d, 
is SOOnm. The incident light is assumed to be a p-wave, 
and the component of the wave vector parallel with the 
film surface is in the x direction, which is parallel to one 
of the lattice vectors. 

Figure 4 shows the spectra of (a) zero-order transmis¬ 
sion, (b) zero-order reflection and (c) total absorption. 
These results were obtained by numerical simulation 
based on the RCWA method. The panel (d) shows the 
dispersion relation of the quasi-bound modes (red dash- 
dotted and blue dashed lines) together with the Rayleigh 
anomaly (dotted lines) and the dispersions of SPPs es- 
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Figure 4. (Color online) Zero-order transmission spectra (a), 
zero-order reflection spectra (b), and total absorption spectra 
(c) for a triangular lattice of nanohole array. Gray dashed 
and dash-dotted lines indicate the dispersion of bound surface 
plasmons (BSPs) and gray dotted lines indicate the condition 
of the Rayleigh anomaly. The bonding (anti-bonding) BSPs 
are shown by the red dash-dotted (blue dashed) lines in (d) 
together with the black solid lines for the empty-lattice bands. 


timated using empty-lattice approximation (black solid 
lines). The dispersion of the quasi-bound modes and the 
Rayleigh anomaly are also shown in the panels (a)-(c) by 
gray dashed (or dash-dotted) lines and gray dotted lines 
respectively. At this point, we can clearly see the coinci¬ 
dence between resonant peaks or dips in the spectrum 
and the dispersions of quasi-bound modes. Although 
there is an abrupt change on the line of the Rayleigh 
anomaly, the profile of this change is still quite different. 
Hence, we conclude that the resonant structure of spectra 
is produced by quasi-bound modes. 

Furthermore, we can also see that there are crossing 
points between red and blue lines. This means that there 
are two kinds of modes which cannot couple with each 
other due to the difference in symmetry. Since these two 
modes contribute independently to the absorption, it is 
expected that the total absorption would become high at 
the crossing points. Indeed, at the crossing point between 
the blue line passing through the point F and the red line 
passing though the point B, the total absorption rate 
becomes quite high, to about more than 90% as shown 
in the panel (c). 
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Figure 5. (Color online) Electric field distributions in the xz 
plane for an anti-bonding mode (a), and a bonding mode (b). 


B. Bonding and Anti-bonding Bound Modes 

As we pointed out in the last subsection, there are 
two types of branches as represented by the red dash- 
dotted lines passing through points A, B and C and blue 
dashed lines passing through points D, E and F in the 
panel (d) of Figure 4. The blue dashed lines are almost 
in accordance with the black solid lines, which indicate 
the empty-lattice bands of SPPs, while the slope of the 
red dash-dotted lines are rather shallow. 

Figure 5 presents the electric field distributions in the 
xz plane (a) at point F and (b) at point A as shown 
in the panel (d) of Figure 4. From the panel (a), we 
see that the charge distribution is anti-symmetric about 
the z = 0 plane, and the intensity of the electric field 
strengthens mainly outside of the metal surfaces between 
the holes. On the other hand, from the panel (b), we 
see that the charge distribution is symmetric about the 
z = 0 plane and the intensity of the electric field concen¬ 
trates around the holes, specifically inside of the metal 
film region. These features indicate that the modes on 
the blue lines are anti-bonding (AB) modes in which the 
SPPs of both sides of the gold film are combined anti¬ 
symmetrically like a long-range SPP of a thin metal film 
[27]. On the other hand, the modes on the red lines 
are bonding (B) modes like a short-range SPP of a thin 
metal film. The symmetry of each bound mode can be 
checked clearly using the CM method. The dispersion re¬ 
lations of B-SPP and AB-SPP calculated using the CM 
method (Eq. (2)) under the symmetric {A = B) and anti¬ 
symmetric (A = —B) conditions are shown in Figures 7 
and 8 , respectively. 

In the AB-SP where the electric field is expelled from 
the hole area, each SPP propagates on each surface with¬ 
out feeling much of the holes. This in turn makes the dis¬ 
persions sit near the empty-lattice bands. We can also 
see that the resonant structure of the spectra created by 
AB-SPP is narrow due to the darkness of the mode. On 
the other hand, the B-SPP has large contribution from 
the waveguide modes in the holes and the dispersions are 
shifted largely from the empty-lattice bands. The shal¬ 


lowness of the band indicates that the B-SPP is much 
like a cavity mode whose electric field is confined in the 
nanoholes. 

The lowest energy B-SPP mode can also be consid¬ 
ered as a coupled spoof surface plasmon, since the disper¬ 
sion saturates around the cutoff frequency of the HEn 
mode, 355 THz. But then what are higher energy B- 
SPP modes? These modes may be considered as coupled 
multipole spoof surface plasmons which are geometrically 
induced coupled surface modes created by the coupling 
between the SPP on the metal surface and the higher 
order waveguide modes of the nanoholes as discussed in 
the next subsection. 


C. Multipole Nature of Bound Modes 

Figure 6 shows the electric field distributions at 20nm 
above the surface of the gold film for the quasi-bound 
modes at the points shown in Figure 4(d). We can 
see monopole (D), dipole (A, E), quadrupole (B, F) 
or hexapole (C) texture for each branch of the quasi¬ 
bound modes. The most striking feature is that a strong 
coupling between adjacent holes seems to exist in the 
quadrupole texture (B, F) as well as positive charges are 
aligned in the a:-direction to form a stripe texture. 

In order to see the role played by the waveguide modes 
to create these multipole textures, we have analyzed the 
contribution of waveguide modes for each branch using 
the CM method. Figure 7 shows the rate of contribu¬ 
tion of the waveguide mode whose name is written at 
the bottom left corner of each panel. The bright color 
means that the rate of contribution of the mode is high. 
It is clearly seen that each branch of the B-SPPs namely, 
A, B and C is mainly created by a single waveguide 
mode, HEii, HE 21 or HE 31 respectively. This means 
that the coupled spoof surface plasmon yields higher en¬ 
ergy branches with multipole textures according to the 
higher order waveguide modes. 

The contribution of the waveguide modes to the AB- 
SPPs is more complicated. There is non-negligible con¬ 
tribution from the TM-like modes, i.e. TMqi and EHn. 
Furthermore, a kind of band inversion seems to occur 
between two bands produced by HEn and HE 21 modes. 
The contribution rates of HEn and HE 21 modes are re¬ 
versed around the points E and F. This could be ex¬ 
plained as a result of band anti-crossing arising from the 
band inversion between HEn and HE 21 bands near the 
P point. 

In order to reveal the origin of these band structures, 
we approximately analyze the dispersion relation of the 
quasi-bound mode produced by a single waveguide mode 
either from HEn or HE 21 . This is done using a meta¬ 
material treatment as in [ 2 ], where only a single mode 
a{= HEii or HE 21 ) inside the hole is taken into account 
and only the p-polarized zero-order diffraction mode is 
considered in the quantity Gaa governing the coupling 
between holes. Starting with the mean-field approxi- 
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Figure 6. (Color online) Electric field distributions at 20nm 
above the surface of the gold film for the points shown in Fig. 
4(d). 



Figure 7. (Color online) Rate of contribution of each waveg¬ 
uide mode to create bonding bound modes calculated using 
the coupled-mode method. 
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Figure 8. (Color online) Rate of contribution of each waveg¬ 
uide mode to create anti-bonding bound modes calculated 
using the coupled-mode method. 


mated CM equations (3), we obtain the condition for the 
existence of B-SPP mode as {paa + = 0 


where S = 


fcp), Sa-G^ = 


Ya/ {ZsYa — 7 ) and 7 = i/tan ( 90 , 2 / 1 / 2 ). The definitions 
of the symbols can be found in the Appendix. Neglect¬ 
ing the imaginary part of the permittivity of gold Era, the 
dispersion relation of the bonding quasi-bound mode is 
given by 




(1-|5P) + 


J ■ 


(4) 

From Eq. (4), we can see that the dispersion relation is 
reduced to that of the SPP on the fiat gold surface within 


the SIBC approach, kx ee y/ekujs/l - e/em — 


if the coupling between the waveguide mode and the 
diffraction mode is negligibly small (l^p ~ 0). On the 
other hand, as the frequency approaches the cutoff fre¬ 
quency from below, ly/Pal diverges for HE modes since 
HE modes approach TE modes and the effective admit¬ 
tance Ya goes to zero as qaz goes to zero. Correspond¬ 
ingly, 7 goes to 00 in this limit. Although the absorption 
of gold prevents from diverging, this divergence property 
remains and the dispersion relation flattens around the 
cutoff frequency as long as the film remains thin. 

The right panel of Figure 9 shows the dispersion re¬ 
lation of Eq. (4) considering the diffraction effects in 
the empty-lattice limit. We can see that the dispersion 
curves, dotted lines and dashed lines, converge to the cut¬ 
off frequencies of HEn mode (355THz) and HE 21 mode 
(516THz), respectively. In this panel, the empty-lattice 
bands of the SPPs on a flat gold surface are also shown by 
black solid lines. The set of numbers (m, n) indicates the 
diffraction order which mainly contributes to the branch 
pointed by the arrow. The parallel momentum k for the 
diffraction order (m, n) is described by the expression 
k = {kx, 0 ) -|- mbi + 7162 where bi and 62 are the recipro¬ 
cal lattice vectors for the triangular lattice of the period 

d given by 61 = ^ ^ (o, . There¬ 

fore, the branches for the diffraction orders ( 1 , 0 ), ( 1 , 1 ), 
(—1,0) and (-1,-1) are composed of the EM wave with 
a period equal to pSd, such that 


= 2e‘('=-±^)"cos(^^y^ . 


Also, the branch for the diffraction orders (0,-1) and 
(0,1) is composed of the EM wave with a period equiva¬ 
lent to p3dj2, where 


= 2e'^^-cos(^y 

\V3d 

These periodicities in the y-direction are considered to 
play a major role to produce the bands of B-SPPs. The 
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Figure 9. (Color online) Dispersion relations of the bond¬ 
ing quasi-bound modes produced by each single waveguide 
mode. The left panel shows those of HEn (magenta circles) 
and HE 21 (cyan triangles) modes. The black solid lines and 
the green dash-dotted lines indicate those of 6 modes and 
two modes. The right panel shows those of HEn (dotted 
line) and HE 21 (dashed line) modes obtained approximately 
with metamaterial treatments. The black solid lines indicate 
empty-lattice bands of the SPP on a flat gold surface. The 
magenta (cyan) solid lines indicate the branches selected to 
form the quasi-bound modes of HEn (HE 21 ) as a guide for 
the eye. 


left panel of Figure 9 shows the dispersion relations calcu¬ 
lated by the mean-field approximated CM equations (3) 
including only a single waveguide mode either HEn (ma¬ 
genta circles) or HE 21 (cyan triangles) but with the pres¬ 
ence of as many diffraction modes as needed to achieve 
convergence. Taking into account the fact that the dis¬ 
persion relation moves closer to the light line if higher- 
order diffraction modes are included [28, 29], it is reason¬ 
able to think that the branches of B-SPP bands created 
by HEn mode are composed of two types of branches: 
(i) the \f^d periodic branches from the flat-surface SPP 
bands that correspond to the black solid lines in the 
right panel and (ii) the •\/3d/2 periodic branches from the 
bands of coupled spoof surface plasmons that represent 
the dotted lines in the right panel. The diffracted waves 
with \/id periodicity may couple strongly with the dipole 
lattices created by HEn mode. In effect, coupled spoof 
surface plasmons are formed and dispersion relation are 
put downward since positive and negative charges align 
alternately in the y-direction to form \/id periodic tex¬ 
ture as you can see in Figure 6 A. Besides, the diffracted 
waves with \fZdl2 periodicity may couple strongly with 
the quadrupole lattices created by HE 21 mode since pos¬ 
itive or negative charges produce bond between adjacent 
holes to form stripe texture with •\/3d/2 periodicity as 
you can see in Figure 6 B. Considering the couplings be¬ 
tween HEii and HE 21 or with higher order waveguide 
modes, the texture seems to be optimized and yield sim¬ 
pler band structure (see green dash-dotted lines and black 
solid lines in the left panel). 

The dispersion relation of the AB-SPP mode is ob¬ 
tained by replacing 7 with 7 “^ in Eq. (4) as shown in the 
right panel of Figure 10. In this case, the influence of the 
coupling between waveguide mode and diffraction modes 



Figure 10. (Color online) Dispersion relations of the anti¬ 
bonding quasi-bound modes produced by a single waveguide 
mode. The meaning of each line is the same as in Figure 9. 
The lines at the right panel correspond to a thicker film with 
thickness of 150nm. 



Figure 11. (Color online) Fano structures in the transmission 
(left panel) and reflection (right panel) spectra. 


become weak as h becomes small. The results for the film 
with h = 150nm are shown in figure for clarity. The left 
panel shows the dispersion relations of AB-SPP obtained 
with the same procedure as used in the left panel of Fig¬ 
ure 9. These results are seemingly produced by the same 
mechanism as those of B-SPP: only the 'JZd {^JidjT) pe¬ 
riodic branches are put downward due to the periodicity 
of the dipole (quadrupole) texture for the HEn (HE 21 ) 
mode. Therefore, the multipole textures are the origin of 
the band inversion and anti-crossing between HEn and 
HE 21 bands. 


D. Fano Resonance 

As mentioned above, the AB-SPPs and the higher en¬ 
ergy B-SPPs have multipole nature, which can be con¬ 
sidered as dark modes. The peak and dip structures in 
the spectra would be attributed to the Fano resonances 
between the brighter modes and the darker modes. In¬ 
deed, if the loss of gold is low to about 10 % of the normal 
loss, then we can clearly see a typical asymmetric Fano 
structures that consist of pairs of peak and dip at nearby 
frequencies [14, 15] as shown in Figure 11. 

Thus, the complicated and distinct structure in the 
transmission and reflection spectra of this system can be 
attributed to the Fano resonances between various orders 
of multipole surface plasmons. 

























































IV. CONCLUSION 

By thorough numerical simulations using RCWA and 
CM methods, we have found that the higher order 
waveguide modes of nanoholes combine the two SPPs 
at both sides of metal him to create two types of quasi¬ 
bound EM modes. These are the (i) AB-SPP with anti¬ 
symmetric charge distribution and (ii) B-SPP with sym¬ 
metric charge distribution about the center plane of the 
him. Due to the multipole nature of the waveguide 
modes, the electric held distribution near the him surface 
forms a multipole texture. The Fano resonances between 
various orders of these multipole surface plasmons pro¬ 
duce sharp peak-dip structure in the transmission and 
rehection spectra. 

The B-SPP modes can be considered as coupled mul¬ 
tipole spoof surface plasmons whose electric helds are 
mainly conhned in the nanohole region and the waveguide 
modes play a major role in determining their dispersions. 
In other words, the nanoholes act like artihcial atoms and 
the multipolar EM wave in a hole hops to other hole via 
SPPs just like the electron in a lattice of real atoms that 
have higher atomic orbitals. More so, the dispersion of 
the AB-SPP mode is determined by the combination of 
multiple waveguide modes which yields a hne band struc¬ 
ture due to the band inversion and anti-crossing induced 
by the multipole textures. 

Further study of these two types of multipole surface 
plasmons will open the possibility to create various artifi¬ 
cial materials with high functionality utilizing their high 
degrees of freedom, e.g. a perfect absorber using multiple 
crossing points of the band structure or a high-resolution 
bio-sensor using multiple resonant transmissions. 


Appendix A: Coupled-Mode Method for Real Metal 


piecewise constants whose values in the tth homogeneous 
medium are given by 


w(d _ 

-'aTE “ 


1 Qctz 
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(A4) 


where Zq and are the impedance and the wavenumber 
in the vacuum, and ni are the relative permittivity and 
relative permeability in the ith homogeneous medium, 
and Qaz is the propagation constant of the mode a. This 
representation indicates that the magnetic field of a mode 
is constructed by a position-dependent admittance from 
the corresponding electric field. If we use Dirac’s notation 
to describe the electric field, such that 


= (f^a), (A5) 

the effect of the position-dependent admittance is ex¬ 
pressed by the admittance operator Y which maps the 
electric field to the corresponding magnetic held, such 
that 


-Bz X Ha{r) = (r 
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The orthogonality condition for the modes can be ex¬ 
pressed by the admittance operator as 
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In general, the mode helds of a waveguide that is 
composed of a piecewise homogeneous media can be ex¬ 
pressed as a linear combination of TE and TM compo¬ 
nents, such that 

Ea if) = aa {r)EaTE (r) + ba {r)EaTM if ), (Af) 
-Bz X Ha{r) = aair)YaTE{r^EaTE{r) 

+ ba{r)YaTMir)EaTMir), (A2) 

where Ea and Ha are the electric and magnetic held 
vectors that are respectively projected onto the xy-plane. 
Also, Bz is a unit vector along the z-axis and r = (a;,y) is 
the position vector in the xy-plane. Here, the mode index 
a represents the full information of the waveguide mode 
of a nanohole in the metal him region, such as the “HEn 
horizontal mode” [5]. It may also represent the parallel 
momentum k and the polarization a{= p or s) for plane- 
wave modes in the homogeneous dielectric layers where 

(Hs^hs) = (1>0) or (^fep^^p) = (0>1)- The position- 
dependent admittances for the TE and TM modes are 


where 6ap is the Kronecker delta and * denotes the com¬ 
plex conjugate. 

Using this dehnition of the admittance operator, the 
CM equation (2) can be derived in a similar manner as 
the original derivation [19]. The a component of the 
coefHcient vector I and the a/3 component of the matrix 
that appeared in Eq. 2 are expressed as 



(A9) 

(AlO) 


Here, we introduced an operator = 1 ± ZgY which 
has an expectation value = 1 ± Z^Y^^ for plane waves 
using the surface impedance 
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with and /Xm being the relative permittivity and the 
relative permeability of the metal. Additionally, the 
transmission and reflection coefficients are expressed as 


(A12) 
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k(7 r-\- ^kko ^^0 

ka 


+ ^^{(ka / ajAa + (ka /+ , 


kcr a 


(A13) 


where Aa (Ba) is the coefficient for the mode a prop¬ 
agating in the +z {—z) direction at the top (bottom) 
interface, and h denotes the him thickness. 

If we adopt a kind of mean-held approximation for the 
admittance operator, such that 


at the top and bottom interfaces as 

(A16) 

A(, = A„e'«-V(r + i3o/a, (A17) 

then the CM equations can be expressed as Eq. (3). The 
a/3 component of the matrices G, S and are given 
as 


< \ ^ ^k<7 j 
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The details of the derivation will be published in an¬ 
other avenue. 
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and introduce the modal amplitudes of the electric held 


ACKNOWLEDGMENTS 

This work was supported in part by KAKENHf No. 
23651149 from MEXT of Japan. 


[1] T. W. Ebbesen, H. J. Lezec, H. F. Ghaemi, T. Thio, and 
P. A. Wolff, Nature 391, 667 (1998). 

[2] F. J. Garcia-Vidal, L. Martm-Moreno, T. W. Ebbesen, 
and L. Kuipers, Rev. Mod. Phys. 82, 729 (2010). 

[3] P. Lalanne, J. G. Rodier, and J. P. Hugonin, Journal of 
Optics A: Pure and Applied Optics 7, 422 (2005). 

[4] H. A. Bethe, Phys. Rev. 66, 163 (1944). 

[5] A. Roberts, J. Opt. Soc. Am. A 4, 1970 (1987). 

[6] J. B. Pendry, L. Martm-Moreno, and F. J. Garcia-Vidal, 
Science 305, 847 (2004). 

[7] L. Martm-Moreno, F. J. Garcia-Vidal, H. J. Lezec, K. M. 
Pellerin, T. Thio, J. B. Pendry, and T. W. Ebbesen, 
Phys. Rev. Lett. 86, 1114 (2001). 

[8] Z. Ruan and M. Qiu, Phys. Rev. Lett. 96, 233901 (2006). 

[9] F. J. Garcia-Vidal, E. Moreno, J. A. Porto, and 
L. Martm-Moreno, Phys. Rev. Lett. 95, 103901 (2005). 

[10] R. Wood, Phil. Mag. 4, 396 (1902). 

[11] M. Sarrazin, J.-P. Vigneron, and J.-M. Vigoureux, Phys. 
Rev. B 67, 085415 (2003). 

[12] O. M. L. Rayleigh, Proc. R. Soc. Lond. A. 79, 399 (1907). 

[13] A. Hessel and A. A. Oliner, Appl. Opt. 4, 1275 (1965). 

[14] C. Genet, M. van Exter, and J. Woerdman, Optics Com¬ 
munications 225, 331 (2003). 

[15] B. luk’yanchuk, N. 1. Zheludev, S. A. Maier, N. J. Halas, 
P. Nordlander, H. Giessen, and C. T. Chong, Nature 
Mater. 9, 707 (2010). 

[16] The contributions of multipolar resonances to the trans¬ 
mission suppression effect in a perforated ultrathin 
metallic film have been investigated in M. Liu, Y. Song, 


Y. Zhang, X. Wang, and C. Jin, Plasmonics 7, 397 

( 2012 ). 

[17] L. Li, J. Opt. Soc. Am. A 14, 2758 (1997). 

[18] T. Weiss, Advanced numerical and semi-analytical scat¬ 
tering matrix calculations for modem nano-optics, Ph.D. 
thesis, University of Stuttgart (2011). 

[19] F. de Leon-Perez, G. Brucoli, F. J. Garcia-Vidal, and 
L. Martin-Moreno, New J. Phys. 10, 105017 (2008). 

[20] A. Taflove and S. C. Hagness, Computational Electro¬ 
dynamics: The Finite-Difference Time-Domain Method, 
3rd ed. (Norwood, MA: Artech, 2005). 

[21] A. F. Oskooi, D. Roundy, M. Ibanescu, P. Bermel, J. D. 
Joannopoulos, and S. G. Johnson, Computer Physics 
Communications 181, 687 (2010). 

[22] J. B. Pendry, D. Schurig, and D. R. Smith, Science 312, 
1780 (2006). 

[23] U. Leonhardt and T. Philbin, Geometry and Light: The 
Science of Invisibility (Dover Publications, 2010). 

[24] See Supplemental Material for additional figures. 

[25] J. D. Jackson, Classical Electrodynamics 3rd ed. (Wiley, 
New York, 1999). 

[26] J. R. Taylor, Scattering Theory (John Wiley & Sons, New 
York, 1972). 

[27] D. Sarid, Phys. Rev. Lett. 47, 1927 (1981). 

[28] F. J. Garcia de Abajo and J. J. Saenz, Phys. Rev. Lett. 
95, 233901 (2005). 

[29] E. Hendry, A. P. Hibbins, and J. R. Sambles, Phys. Rev. 
B 78, 235426 (2008). 



